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We provide a rigorous solution to the problem of constructing a structural evolution for a network 
of coupled identical dynamical units that switches between specified topologies without constraints 
on their structure. The evolution of the structure is determined indirectly, from a carefully built 
transformation of the eigenvector matrices of the coupling Laplacians, which are guaranteed to 
change smoothly in time. In turn, this allows to extend the Master Stability Function formalism, 
which can be used to assess the stability of a synchronized state. This approach is independent from 
the particular topologies that the network visits, and is not restructed to commuting structures. 

Also, it does not depend on the time scale of the evolution, which can be faster than, comparable 
to, or even secular with respect to the the dynamics of the units. 

PACS numbers: 89.75.He, 05.45.Xt, 87.18.Sn, 89.75.-k 


Networked structures, in which sets of distributed dy¬ 
namical systems interact over a wiring of links with non¬ 
trivial topology, are a key tool to investigate the emer¬ 
gence of collective organization in many systems of inter¬ 
est EH. The analysis of synchronized states is particu¬ 
larly relevant, as they play a crucial role in many natural 
systems, such as brain networks Q, or ecological com¬ 
munities Q. In the past decade, the emergence of syn¬ 
chronized states has been extensively reported and stud¬ 
ied with a notable emphasis on the effect of complex 
static topologies on synchronization properties (7Hl7l|. 
Nonetheless, static network models do not adequately 
describe the processes that arise because of mutations in 
some biological systems, such as infectious bacterial pop¬ 
ulation, which are known to have adaptive mutation rates 
that can become very highly elevated El El- These 
require, instead, the use of time-dependent topologies 
whose evolution occurs over time scales that are commen¬ 
surate with those of the node dynamics [ 23 , 1 m. ^ power¬ 
ful tool to assess the stability of the synchronous solution 
in networks of N identical nodes with diffusive coupling, 
is the so-called Master Stability Function (MSF) Q. If 
the evolution of such networks is along structures whose 
Laplacians commute at all times, synchronization can be 
significantly enhanced. In fact, studies have shown that 
it can be achieved even when the connection structure 
at each time would not allow a synchronized state, in 
the static case [l^, 1131 . The far more realistic situation 
of networks whose coupling matrices do not necessarily 
commute has also generated significant interest. This has 
led to the development of several results about synchro- 
nizability in systems with specific properties, including 
the study of the synchronous state in the case of moving 
neighbourhood networks [mi: the rigorous derivation of 
sufficient conditions for synchronization in fast switching 


networks ^-nd the analysis of the system dynamics 
in the so-called blinking limit 1^. In this Article, 
we propose a framework that allows to assess the stabil¬ 
ity of a synchronous solution by means of an MSF when 
the evolution of the topology is fully general and uncon¬ 
strained, and without assuming any hypothesis on the 
time scales of the topological evolution. 

To this purpose, let us consider a network of N iden¬ 
tical systems, evolving according to 

N 

Xi = f (x,) (t) h (Xj) . (1) 

t=i 

Here dot denotes time derivative, x^ is an m-dimensional 
row vector describing the state of the node, a is 
the interaction strength, and f and h are two vecto¬ 
rial functions describing the local dynamics of the nodes 
and the output from a node to another, respectively. 
Also, L{t) = S (t) — W (t) is the Laplacian of the net¬ 
work describing the time evolution of the connections. 
In this expression, W is the weighted adjacency matrix, 
and S is the diagonal strength matrix of the network: 
Sij = Sij To fix the ideas, assume that the 

network has an initial structure, with Laplacian Lq, that 
is constant from time t = 0 to t = to. Also assume 
that, at time ti > to, the network is found in a con¬ 
figuration with Laplacian Li ^ Lq. Note that in the 
following we consider Laplacians with Hermitian struc¬ 
ture. Also, as L (t) is a Laplacian matrix, the sum of 
each row vanishes, its diagonal elements are strictly posi¬ 
tive and its off-diagonal elements are non-positive. Now, 
let Xg be the state vector indicating the synchronized 
solution, and define the miV-dimensional column vector 
i5X = (5xi,... , 6xn)^ , representing the global deviation 
from the synchronized state. From Eq. [TJ to linear order 
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in SX, one has 

= (1 (g) Jf (x,) - ctL (t) (g) Jh (x,)) SX , (2) 

where 1 is the A^-dimensional identity matrix, (g) denotes 
the direct product, and J is the Jacobian operator. The 
vector SX can be written at each time as a sum of di¬ 
rect products of the eigenvectors of L (t) and a time- 
dependent set of N ?7i-dimensional row-vectors rj^ (t): 

N 

JX = ^ Vi (t) (g) J7, (t) . 

Multiplying vj from the left to both sides of Eq. [5J one 
gets 


_d 

dt 





( 3 ) 


where, for the sake of brevity, we defined Kj = 
(Jf (Xg) — Jh(xs)) and i^j = aXj (t), in which Xj (t) is 
the eigenvalue of L (t). To assess the stability of this 
dynamical system, one can compute the Master Stability 
Function, which represents the dependence on v of the 
largest Lyapunov exponent Amax associated to the equa¬ 
tions [3J Then, the stability criterion for a given v is that 
the time averages of A^ax in the direction of all eigen¬ 
vectors are negative. This allows to study systems with 
highly non-trivial behaviour. As an example, one can 
consider an evolving system where each “frozen” connec¬ 
tion topology has at least one direction in which Amax is 
positive, but the synchronization manifold is transversely 
stable [2^. Similarly, one can detect instabilities intro¬ 
duced by the evolving nature of the Laplacian in systems 
where the synchronization manifold in each frozen config¬ 
uration is attracting [2^. This is particularly useful, as it 
is well-known that non-linear perturbations in a system 
can destroy the stability of the synchronized state [2^. 

Note that using the MSF is not the only possible 
method to assess the stability of a synchronized state. 
For instance, one could construct the Lyapunov function 
for the synchronization manifold. This would guarantee 
the stability since it is a necessary and sufficient condi¬ 
tion, while in the general case the MSF is only a necessary 
one. However, while it is certainly possible to build the 
Lyapunov function in some specific cases [^ . a general 
construction method is not known. Also, it should be 
noted that in the absence of fixed points or other attract¬ 
ing sets far from the synchronization manifold, the MSF 
provides a sufficient stability condition as well, thereby 
becoming a widely used approach . 

To use the MSF method, one must first note that the 
boxed term in Eqs. [3] explicitly depends on the time vari¬ 
ation of the eigenvectors of L. If the Laplacians Lq and 
Li commute, one can choose to study the problem in the 
common basis of eigenvectors. In this reference frame 
the eigenvectors do not change. Thus, the boxed term 


vanishes and Eqs. |3] reduce to a set of variational equa¬ 
tions. However, if we allow the network to switch to a 
different structure without imposing the extra require¬ 
ment of commutativity, the eigenvector variation must 
be taken into account. Note that this forbids instan¬ 
taneous jumps between the two structures, as a sudden 
change in the eigenvectors would cause their derivatives 
to be not defined. Therefore, the goal is constructing a 
smooth evolution process from t = to to t = ti to allow 
the system to evolve between between the two topolo¬ 
gies while keeping the eigenvector elements differentiable. 
To achieve this, we first consider the matrices A and B, 
consisting of the eigenvectors of Lq and Li, respectively, 
and describe how to transform one into the other via a 
proper rotation around a fixed axis. Then, we use this 
framework to find a transformation between Lq and Li. 
The rotation evolving A into B takes the form of a one- 
parameter transformation group Gg such that GqA = A 
and GiA = B. Note that, as we will use this mapping 
to build the Laplacian, we must also impose the extra 
requirement that the vector a, corresponding to the null 
eigenvalue, is kept constant by Gg for all 0 ^ s ^ 1. 

In general, the transformation O from A to H is a 
rotation, which can be found by solving the linear system 
of equations in N'^ variables OA = B. It is convenient 
to work in the basis defined by the matrix A. In this 
basis, A = 1. Without loss of generality, assume that 
the conserved vector a is the first vector of A. Then, the 
transformation matrix O has the form 



/I 

0 

0 • 

■ \ 


0 

Oil 

023 ■ 

• O 1 N 

o = 

0 

032 

033 ■ 

■ OoN 


VO 

OAr 2 

Oats ■ 

■ OnnJ 


As O G O (N), it is a proper rotation if its determinant 
is 1, or a composition of rotations and reflections if its 
determinant is —1. The determinant of O equals the 
determinant of the minor O' obtained by removing the 
first row and the first column from O. Thus, we only need 
to find a solution to the problem in — 1 dimensions. 
We will henceforth use primes when referring to objects 
in A^ — 1 dimensions. 

From the considerations above, it is Gq = V G'^ = O', 
and, of course. O' G 0{N — \). Thus, the problem 
is equivalent to determining the possibility of finding a 
path between the identity and O'. If |0'| = 1, then 
O' G SO [N — 1). But SO [N — 1) is the connected iden¬ 
tity component of the orthogonal group, and additionally, 
since it is a manifold, it is path-connected (j^ . Thus, for 
every orthogonal {N — 1) x {N — 1) matrix in it, there is a 
continuum of orthogonal matrices of the same dimension 
connecting it to the identity. Each point along this path 
corresponds to an orthogonal matrix that can be embed¬ 
ded in SO {N) by adding a 1 in the top left corner. Since 
every such embedded matrix keeps the synchronization 
manifold vector invariant, a parametrization of the path 
provides a direct solution to the original problem. 
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If, instead, \0'\ = —1, then O' S 0{N — 1) \ 
SO [N - 1). While O [N - l)\SO {N - 1) is also a con¬ 
nected topological space, the identity does not belong 
to it [ 3 ^. Thns, no path connects the identity to O'. 
However, in our case, the labeling of the vectors is irrel¬ 
evant. In other words, provided that the vector a is kept 
constant, one can arbitrarily swap two vectors in the ba¬ 
sis given by the matrix A, obtaining a new matrix C. Of 
course, this imposes a swap of the corresponding columns 
in the transformation matrix O as well. But swapping 
two rows in a matrix changes the sign of its determinant. 
This means that the new matrix O' S SO {N — I), and a 
path connecting it to the identity can be found. The only 
consequence of the swap is a change of the order of the 
eigenvalues. As we are considering unlabelled networks, 
the problem can always be solved. To build an explicit 
solution, we factor the transformation O' into rotations 
and reflections in mutually orthogonal subspaces. 

Before describing the actual procedure, we recall a use¬ 
ful, twofold result. First, any orthogonal operator X in 
a normed space over R induces a I- or 2-dimensional in¬ 
variant subspace. To find one such subspace, first define 
the unitary operator U acting on a -f ib as ?7 (a -|- ib) = 
Xa + iXb, where a and b belong to Then, find a 
non-vanishing eigenvector x = xr -|- ixi of U. The span 
of xr and xi defines the invariant subspace: applying X 
to any linear combination of xr and xj produces a vector 
that is still a linear combination of xr and xi. Also, if 
the corresponding eigenvalue is complex, xr and xj are 
orthogonal. 

Using this, we can describe the following algorithmic 
procedure to build a transformation O of A into B: 

1. Express the problem in the basis A, in which O has 
the form of Eq. SI 

2. Consider the operator O' obtained from O by re¬ 
moving the first row and the first column, and let 
d be its dimension. 

3. Build the operator U that acts on a -|- ib as 
U {a + ib) = O'a + iO'b, where a and b belong 
to 

4. Find an eigenvector x = xr -|- ixi of U, with eigen¬ 
value A. 

5. Normalize xr and xj. 

6. If A £ R then 

(a) Pick the non-vanishing component between 
Xr and xj. If both are non-zero, choose one 
randomly. Without loss of generality, assume 
this is Xr. 

(b) Create d—1 other orthonormal vectors, all or¬ 
thogonal to Xr, and arrange all these vectors 
so that Xr is the last of them. This set of 
vectors is an orthonormal basis C of R^^. 


(c) Change the basis of the d-dimensional sub¬ 
problem to C. In this basis, all the elements 
in the last row and in the last column of O' 
will be 0, except the last one, which will be 
± 1 . 

(d) If d > 1, consider a new operator O' obtained 
from the old O' by removing the last row and 
the last column. Let d be its dimension, and 
restart from step 3. Otherwise stop. 

If instead A ^ R, then 

(a) Create d — 2 other orthonormal vectors, all 
orthogonal to xr and xi, and arrange all these 
vectors so that xr and xj are the first two of 
them. This set of vectors is an orthonormal 
basis C of R'^. 

(b) Change the basis of the d-dimensional sub¬ 
problem to C. In this basis, all the elements in 
the first two rows and in the first two columns 
of O' will be 0, except the first two. 

(c) If d > 2, consider a new operator O' obtained 
from the old O' by removing the first two rows 
and the first two columns. Let d be its dimen¬ 
sion, and restart from step 3. Otherwise stop. 

At each iteration of the steps above, I or 2 dimen¬ 
sions are eliminated from the problem. All the subse¬ 
quent changes of base leave the already determined ele¬ 
ments of O unchanged, because they act on orthogonal 
subspaces to those already eliminated. The procedure re¬ 
constructs O piece by piece with a block-diagonal form, 
in which the blocks correspond to the action of the or¬ 
thogonal operator on the invariant subspaces. If the sub¬ 
space is 1-dimensional, then the block is a single ±1 ele¬ 
ment. If instead the subspace is 2-dimensional, then its 
block is either a rotation or a reflection, i.e., it is either 

/cos a -smaA /^±1 0 Y the 1-dimensional 

l^sina cos a J 0 =Fl/ 

invariant subspaces induced by the operator correspond 
either to leaving one direction untouched, or to reflect¬ 
ing the system about that direction. Conversely, the 2- 
dimensional invariant subspaces correspond to rotations 
in mutually orthogonal planes. 

Once this form of O is found, permute the basis vectors 
from the second onwards so they correspond, in order, 
first to all the actual rotation blocks, then to the —1 ele¬ 
ments, and finally to the -1-1 elements. The new form of 
the transformation matrix. Oat, is simply On = TOT~^, 
where T is the required change-of-basis matrix. 

Note that the determinant oi On could still be — 1. 
However, as seen before, in this case, one can relabel 
two vectors of the original basis, inducing a swap of the 
corresponding columns of On ■ To perform this, first note 
that, if I Oat I = —1, then the number of —1 elements in 
On must be odd. Then, there are three possible cases. 

If 0 )y has at least one -1-1 element, swap the basis 
vectors corresponding to the first —1 and the first -|-1 
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elements in O'. 


cos” blocks is 


Then, the first block after the “sin- 
^ ' Swapping the labels of the the 


1 

0 -1 


two corresponding vectors, makes the block in O' become 


0 

-1 


_/cos(-f) -sin(-|) 




sm(-f) cos (-- 2 ; 

If instead has no +1 elements and only one —1 el¬ 
ement, the basis vectors to swap are those corresponding 
to the first two vectors in the last 3x3 block of O', that 
becomes 


( —sindfe cosdfe 0\ 

cosr^fc sindfe 0 I 
0 0 -ij 


hrst two after the “sin-cos” blocks. Their 3x3 block is 
now 


M = 




Next, do a change of basis as described in the previous 
case. The eigenvector matrix of M is 


V = 




Now, perform one more basis change, leaving all the basis 
vectors unchanged, but mapping the last three into the 
eigenvectors of M, which are 

0 (—sec —tan sin ^ (sec — tan ) sin^ -t-\ 

0 1 1 / 
10 0 / 

The new form of M becomes 

/-I 0 0\ 

M' = V'^MV = 0 -10 . 

\o 0 ij 

Finally, if has no -|-1 elements and at least three 
— 1 elements, swap the basis vectors corresponding to the 

_I 


so, after the basis change, the new form of M is once 
more 

/-I 0 0\ 

M' = V'^MV = 0 -10 . 

\o 0 ij 

Regardless of the original value of IOat], take now ev¬ 
ery two subsequent —1 elements, if there are any, and 

change their diagonal block into sin7r\ cpjjjg 

° ° f^SlUTT COSTT J 

yields the hnal general form for the transformation ma¬ 
trix, which can be turned into the required transforma¬ 
tion group via the introduction of a parameter s G [0,1]: 


/I 

0 

0 


0 

0 

0 

0 

0 

0 

0 

•• 0 0 \ 


cos('j9is) 

— sin('i9i 5 ) 


0 

0 

0 

0 

0 

0 

0 

•• 0 0 \ 

0 

sin('j9is) 

cos('i9i s) 


0 

0 

0 

0 

0 

0 

0 

■■00 





0 

0 

0 

0 

0 

0 

0 

■■00 

0 

0 

0 

0 

COs(l9fc5) 

— sin(i9fcs) 

0 

0 

0 

0 

0 

••00 

0 

0 

0 

0 

sin(i9fcs) 

COs(l?fe5) 

0 

0 

0 

0 

0 

■■00 

0 

0 

0 

0 

0 

0 

cos(-f s) 

-sin(-f s) 

0 

0 

0 

■■00 

0 

0 

0 

0 

0 

0 

sin(-f s) 

cos(-f s) 

0 

0 

0 

-- 0 0 

0 

0 

0 

0 

0 

0 

0 

0 

cos(7rs) 

— sin(7rs 

) 0 

■■00 

0 

0 

0 

0 

0 

0 

0 

0 

sin(7rs) 

cos(7rs^ 

0 

■•0 0 


6 

6 

6 

6 

6 

6 
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6 

6 

6 

■' ib j 

^0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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Notice that when s = 0, Gs = Go = 1, and when 
s = 1, Gs = Gi = O. Also, for every value of s, the 
determinant of Gs is always 1, and the first vector is 
kept constant, which means that Gg describes a proper 
rotation around the axis defined by the eigenvector cor¬ 
responding to the synchronization manifold. Moreover, 
the first vector has always been left untouched by all the 
possible basis transformations. Thus, as s is varied con¬ 
tinuously between 0 and 1, the application of Gg sends 
A into B continuously, as needed. 

Finally, to describe how to obtain a transformation 


between the two Laplacians Lq and Li, let R be the 
change-of-basis matrix resulting from the composition of 
all basis changes done in step 6 of the method, the per¬ 
mutation of the vectors to obtain O, and the possible 
hnal adjustment in case of negative determinant. Also, 
to simplify the formalism, in the following we let Bq = A 
and Bi = B. Since Bq and Bi are matrices of eigenvec¬ 
tors, it is BJLqBq = Dq and BJLiBi = Di, where Dq 
and Di are the diagonal matrices of the eigenvalues of 
To and Li. Also, as Gg is a group of orthogonal transfor¬ 
mations, it is GsRBq = RBs VO ^ s ^ 1, where Bg is a 
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basis of . Multiplying this by on the left yields 

R^GsRBo = R^RBs = B,, (6) 

where we used R^ = R~^. 

Now, note that all the basis changes performed are 
between orthonormal bases. Thus R is an isometry, as 
is any Gs, since they are all proper rigid rotations, and 
any Bg defines an orthonormal basis of . Then, the 
Laplacian for the parameter s is given by the matrix Lg 
that solves the equation 

BjLgBg=Dg, (7) 

where Dg is a diagonal matrix whose elements are the 
eigenvalues of Lg, and Bg consists of the eigenvectors 
of Lg. However, the equation above has two unknowns, 
namely Lg and Dg. 

This provides a certain freedom in describing the evo¬ 
lution of the eigenvalues of the Laplacians. For instance, 
one can choose the simplest evolution, which is given by 
a set of linear transformations. Then, for all 1 ^ i ^ iV, 

= il-s)X^^ +sX^K (8) 

where A® and A^^^ are the eigenvalues of Lq and Li, 
respectively. Note that this allows for the possibility of 
degeneracy of some eigenvalues for some particular value 
of the parameter s*. However, the transformation we 
described leaves all the eigenvectors distinct and sepa¬ 
rate throughout the evolution. Thus, the Laplacian can 
always be diagonalized for any value of s. Multiplying 
Eq. [7] on the right by Bj yields BjLgBgBj = DgBj, 
hence BjLg = DgBj. Multiplying this on the left by 
Bg, it is BgBjLg = BgDgBj, hence Lg = BgDgBj. 
Substituting Eq. |6] into this last equation gives 

Lg = R^GgRBoDg {R^GgRBof 
= R^GgRBoDgB'^R^GjR. 

In the equation above, Bq is known, R and Gg have 
been explicitly built, and Dg is completely determined 

I 


by Eq. [51 Thus, Eq. [S] defines the Laplacian for any given 
value of the parameter s. The evolution of the eigenvalues 
is imposed by Eq. El and the evolution of the eigenvectors 
is given by Eq. El 

It is important to stress that the solution given by the 
linear evolution of the eigenvalues is not unique. There 
could be, in principle, many allowed transformations of 
Lq into Li, each characterized by a specific eigenvalue 
evolution. However, the difference between solutions is 
only in the term in Eq. El since the boxed term does 
not depend on the eigenvalues. Thus, the specific choice 
of eigenvalue evolution could change the phenomenology 
of the system studied, but would not modify how the 
switching between general structures affects the stability. 
This is akin to the results presented in Ref. which 
can be considered a special case of the present treatment 
that occurs when all the Laplacians commute. 

With this result, one can finally describe the system 
evolution through unconstrained topologies. From t = 0 
to t = to, the boxed term in Eq. El vanishes. To compute 
its value during the switch, first note that the eigen¬ 
vector at time t is the column of Bg, with s = . 

’ tl—to 

But then, using Eq. El the fc*** element of the i*** eigen¬ 
vector is 


(v.)fe = m,, = (R^GgRBo)^^ = 

N N N 

= ^ ^ ^ Rrk {Gg)^^ Rqx {Bo)xi ■ 

r—1q—1x—1 


Notice that in the equation above the only term that 
depends on time is (Gs)^,^, since R is just a change-of- 
basis matrix, and Bq is the matrix of eigenvectors of Lq 
at time t = 0. Therefore, it is 


_d 

dt 


N N N 

(vO.-EEE Rrk Rqx 

r—1 q—1 x—1 


1 


ti — to ds 


{Gg) 


rq 


This allows a fully explicit expression for the boxed 
term in Eqs.Elthat accounts for the time variation of the 
eigenvectors: 


N , 

2 = 1 


N N N 


N N 

—^EE 

tl — to ^ ^ 

2=1 k—1 q—^ x—1 


EEE Rrk {Gg)j.q Rqx {Bo)xj 


' N N N 

E/ E/ E/ ^rkRqx {Bq)„, 

r—1q—1x—1 


ds 


{Gs) 


rq 




iM 


However, for all practical purposes, one does not need 
to use the expression above directly. In fact, considering 
that most elements of Gg are 0, it is quite simple to com¬ 
pute and store Bg in a symbolic form. Similarly, most of 
the ^ (Gs)^^ terms are 0. In fact, they vanish if r = 1, if 
g = 1, if |r — g| > 1, if r > 26-1-1, and if g > 26-|-1, where 


6 is the number of “sin-cos” blocks in Gg. Also, for all 
other cases ^ {Gg)xq is proportional to a sine or a cosine. 
Then, define Gg to be the matrix whose (rg) element is 

tiEs {Gs)rq\ also, define Bg = R^GgRBQ. Again, Bg 
can be easily computed and stored in a symbolic form. 
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Figure 1: Network topologies. Each node is a chaotic Rossler oscillator. The edges represent weighted diffusive connections. 
The weights are indicated by labels in the vicinity of the corresponding edges. The network switches between the configurations 
in the two panels, remaining in each for some time. 


Then, Eq. [TO] becomes 

N , N N 

- E E (^) ^* = - E E (^0 M ■ 

i—\ i—1 k—1 

( 11 ) 

Once more, the equation above can be computed symbol¬ 
ically, and evaluated at any particular t, when needed. 
Note that, despite the seemingly complex expressions, 
Eqs. [TO| and HD are straightforward to deal with. This is 
due to the fact that Gg and Gg are always represented 
as tridiagonal matrices, and, as mentioned above, many 
of the non-trivial elements of Gg vanish as well. Thus, 
numerical applications of this approach can benefit not 
only from a restricted amount of needed memory, but 
also from sparse matrix methods that result in a small 
computational complexity. 

The treatment we built is valid for every positive, fi¬ 
nite switching time t* = ti — to between configurations. 
As explained above, this time cannot vanish, lest the 
derivatives in Eq. |3| be not defined. Nonetheless, one 
can wonder about the behaviour of a system when the 
switching time becomes very small, although non-zero. 
To this purpose, first note that this time only appears 
as a multiplicative factor in Eq. jTOj Thus, a very small 
t* would have the effect of making the boxed term in 
Eq. |3| much larger than the purely variational term. In 
this regime, the effects on the the stability of the syn¬ 
chronous solution are due mostly, if not exclusively, to 
the switching process. In other words, if the expression 
in Eq. [TO| yields positive results, the synchronized state 
is made more stable in the corresponding direction, and 
vice versa for negative results, regardless of the contribu¬ 
tion coming from the variational term. Note that this is 


in agreement with the finding that blinking networks can 
greatly facilitate synchronization [Til . . Similarly, one 
can consider the opposite limit, namely that of a secular 
switching for which t* becomes very large, while still re¬ 
maining finite. In this case, for large enough switching 
time, the boxed term in Eq. |3| becomes negligible com¬ 
pared to the rest, and the stability is determined entirely 
by the variational term. This case is very similar to that 
of an evolution along commutative structures (^ . In 
fact, in this regime of quasi-static evolution, the structure 
at any given time t is, to first order, equal to the structure 
at time t-|-dt. Thus, it is [Lt,Lt+dt] = £■ Therefore, one 
can treat this case as the commutative one with the addi¬ 
tion of a small perturbation. Note that this perturbation 
does not change the stability of the synchronized state: 
for a positive variational term, instability is maintained, 
and for a negative one, the synchrony remains stable. 
The only uncertainty happens for the critical condition 
corresponding to a vanishing variational term, for which 
the perturbation can have either effect on stability. 

To illustrate the use of our method, we consider the ex¬ 
ample of a weighted network of = 10 chaotic Rossler 
oscillators, switching back and forth between two topolo¬ 
gies (Fig. [J). Letting the state vector x = {x, y, z), each 
of the oscillators obeys the local dynamics 

f (x) = {—y — z,x + 0.165j/, 0.2 z{x — 10)) , (12) 

with the output function 

h(x) = (0,?/,0) . (13) 

The switching times and the time periods for which the 
network remains in each of the two configurations (per¬ 
manence times) are all set to 0.1. We perform two sim¬ 
ulations, one with interaction strength cr = 1 and one 
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Figure 2: (Color online) Estimating the largest Lyapunov 
exponent. The average of the logarithm of the norm of rj 
(Eq. I14II for the example system converges to approximately 
—0.3 when a — 1 (black solid line), and to approximately 
0.022 when ct = 0.1 (red dashed line), indicating that the 
synchronized state is stable in the first case, and unstable in 
the second. Note the logarithmic time scale. 

with cr = 0.1, to illustrate two different cases and the 
sensitivity of our method. To estimate the largest Lya¬ 
punov exponent associated to the system of Eq. |31 we 
compute the time-average of the logarithm of the norm 
of the vector 

»7= (»72.»73>---.’7jv) (14) 

at each integration step. The value to which < log |J7| > 
converges is Amax- The results, in Fig. [2l show that 
for the cr = 1 case the convergence value is approxi¬ 
mately —0.3, indicating that the synchronized state is 
stable. Conversely, when a = 0.1, the estimated Lya¬ 
punov exponent is just positive, with a value of approxi¬ 
mately 0.022, corresponding to an unstable synchronized 
state. To verify this numerical result, we simulated the 
actual network evolution for the two cases according to 
Eq. [2 and with f and h given by Eqs. [T^ and [T51 above. 
Figure [3] shows the time evolution of the global synchro¬ 
nization error 

1 ^ 

X = (1^* - xi\ + \yi - yi\ + \zi - zi\) . 

(15) 

For the ct = 1 case, after a certain transient, the synchro¬ 
nization error decays to 0. When the interaction strength 
is lowered to cr = 0.1, instead, x eventually starts grow¬ 
ing and oscillates wildly, always taking non-zero values. 
These results indicate that the system is indeed able to 
synchronize in the first case, while it never does in the 
second, in agreement with the numerical calculations of 
Amax- Thus, the simulations not only conhrm the valid¬ 
ity of our treatment, but provide an example for which 
the stability of the synchronized state can be changed 
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Figure 3: (Color online) Stability of the synchronized state. 
The time evolution of the global synchronization error x 
(Eq. I15II shows that the system is able to reach synchroniza¬ 
tion when a = 1 (solid black line), while it never synchronizes 
when CT = 0.1 (dashed red line), confirming the numerical re¬ 
sults. Note the logarithmic time scale and the break in the 
vertical axis. 


by the tuning the parameters controlling the topological 
evolution. 

In summary, we demonstrated how to explicitly solve 
the problem of constructing an appropriate time evolu¬ 
tion of a system of networked dynamical units switching 
between different topologies. Our method builds the evo¬ 
lution from a mapping of the eigenvectors of the graph 
Laplacians of the individual structures, and it ensures 
that the elements of the eigenvectors are differentiable at 
each intermediate time. This enables the use of the Mas¬ 
ter Stability Function for network topologies that evolve 
in time in a fully general and unconstrained way. While 
the connection pathway is not unique, different solutions 
only affect the variational part of the linearized system. 
It has to be remarked that our treatment is valid regard¬ 
less of the time scales involved. There is no restriction 
on the permanence times of the network in each configu¬ 
ration, and the only constraint on the switching times is 
that they do not vanish. In addition, our method intro¬ 
duces a numerical advantage, in that one only needs to 
integrate a set of linear equations coupled with a sin¬ 
gle non-linear one, rather than having to deal with a 
system entirely composed of non-linear differential equa¬ 
tions. Also, this approach does not rely on particular 
assumptions concerning the structures visited by the sys¬ 
tems, and contains the regimes of blinking networks and 
commutative evolution as its limiting cases. Thus, our 
results have a natural application in the study of synchro¬ 
nization events in systems for which the temporal scales 
of the topology evolution are comparable with (or even 
secular with respect to) those characterizing the evolu¬ 
tion of the dynamics in each networked unit. This is a 
common occurrence in many real-world systems, such as 
neural networks, where synchronization can become pos- 











sible due to mutations |^ . or financial market, where 
global properties are affected by adaptive social dynam¬ 
ics. 
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